---
title: "conjoint analysis"
author: "Kota Suechika"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
rm(list=ls())

require(tidyverse)
require(cjoint)
library(cregg)
library(scales)
library(patchwork)
library(extrafont)
require(memisc)
require(stringi)
require(ggpubr)
library(tidytext)

dat <- read_csv("data/cj_data_Lebanon1.0.csv")
dat <- na.exclude(dat)

str(dat)
dat[,c(37:48)] <- data.frame(lapply(dat[,c(37:48)],as.factor))
dat[,c(50:52)] <- data.frame(lapply(dat[,c(50:52)],as.factor))

# Ceasefire: unexpected event during survey
dat$RecordedDate <- as.Date(stri_match_first_regex(dat$RecordedDate, "\\d{4}-\\d{2}-\\d{2}"))

dat$ceasefire <- NA
dat$ceasefire[dat$RecordedDate <= as.Date("2025-01-16")] <- "Pre"
dat$ceasefire[as.Date("2025-01-17") <= dat$RecordedDate] <- "Post"
dat$ceasefire <- factor(dat$ceasefire)

# Task
dat$task <- as.factor(dat$task)

# Gender
dat <- dat |> 
  mutate(Gender = case_when(Gender == 1 ~ "Male",
                            Gender == 2 ~ "Female"))
dat$Gender <- as.factor(dat$Gender)

# Age group
dat <- dat |> 
  mutate(Age_group = case_when(D2 == 1 ~ "18-20",
                               D2 == 2 ~ "21-30",
                               D2 == 3 ~ "31-40",
                               D2 == 4 ~ "41-50",
                               D2 == 5 ~ "51-60"))
dat$Age_group <- as.factor(dat$Age_group)

# Education
dat <- dat |>
  mutate(Education = D3)
dat$Education <- as.factor(dat$Education)

# Income
dat <- dat |>
  mutate(Income = D5)
dat$Income <- as.factor(dat$Income)

# Sect
dat <- dat |> 
  mutate(Sect = case_when(D6 == 1 ~ "Sunni",
                          D6 == 2 ~ "Shia",
                          D6 == 3 ~ "Maronites",
                          D6 == 4 ~ "Druze",
                          D6 == 5 ~ "Ismail",
                          D6 == 6 ~ "Armenian",
                          D6 == 7 ~ "Greek",
                          D6 == 8 ~ "Protestant",
                          D6 == 9 ~ "Minorities",
                          D6 == 10 ~ "Multi-sectarian",
                          D6 == 13 ~ "Secular"))
dat$Sect <- as.factor(dat$Sect)

# Islamism
dat <- dat |> 
  mutate(Islamism = case_when(Q1_1 >= 67 ~ "High",
                              Q1_1 <= 66 & Q1_1 >= 34 ~ "Middle",
                              Q1_1 <= 33 ~ "Low"))
dat$Islamism <- as.factor(dat$Islamism)

# Secularism
dat <- dat |> 
  mutate(Secularism = case_when(Q1_2 >= 67 ~ "High",
                                Q1_2 <= 66 & Q1_2 >= 34 ~ "Middle",
                                Q1_2 <= 33 ~ "Low"))
dat$Secularism <- as.factor(dat$Secularism)

# Zaimism
dat <- dat |> 
  mutate(Zaimism = case_when(Q1_3 >= 67 ~ "High",
                               Q1_3 <= 66 & Q1_3 >= 34 ~ "Middle",
                               Q1_3 <= 33 ~ "Low"))
dat$Zaimism <- as.factor(dat$Zaimism)

# Arabism
dat <- dat |> 
  mutate(Arabism = case_when(Q1_4 >= 67 ~ "High",
                             Q1_4 <= 66 & Q1_4 >= 34 ~ "Middle",
                             Q1_4 <= 33 ~ "Low"))
dat$Arabism <- as.factor(dat$Arabism)

# US
dat <- dat |> 
  mutate(US = case_when(Q2_1 >= 4 ~ "Like",
                        Q2_1 == 3 ~ "Neither",
                        Q2_1 <= 2 ~ "Dislike"))
dat$US <- as.factor(dat$US)

# Iran
dat <- dat |> 
  mutate(Iran = case_when(Q2_2 >= 4 ~ "Like",
                          Q2_2 == 3 ~ "Neither",
                          Q2_2 <= 2 ~ "Dislike"))
dat$Iran <- as.factor(dat$Iran)

# Saudi
dat <- dat |> 
  mutate(Saudi = case_when(Q2_3 >= 4 ~ "Like",
                           Q2_3 == 3 ~ "Neither",
                           Q2_3 <= 2 ~ "Dislike"))
dat$Saudi <- as.factor(dat$Saudi)

# Turkey
dat <- dat |> 
  mutate(Turkey = case_when(Q2_4 >= 4 ~ "Like",
                            Q2_4 == 3 ~ "Neither",
                            Q2_4 <= 2 ~ "Dislike"))
dat$Turkey <- as.factor(dat$Turkey)

# Syria
dat <- dat |> 
  mutate(Syria = case_when(Q2_5 >= 4 ~ "Like",
                          Q2_5 == 3 ~ "Neither",
                          Q2_5 <= 2 ~ "Dislike"))
dat$Syria <- as.factor(dat$Syria)

# Israel
dat <- dat |> 
  mutate(Israel = case_when(Q2_6 >= 4 ~ "Like",
                            Q2_6 == 3 ~ "Neither",
                            Q2_6 <= 2 ~ "Dislike"))
dat$Israel <- as.factor(dat$Israel)

# Axis
dat <- dat |> 
  mutate(Axis = case_when(Q2_7 >= 4 ~ "Like",
                          Q2_7 == 3 ~ "Neither",
                          Q2_7 <= 2 ~ "Dislike"))
dat$Axis <- as.factor(dat$Axis)

# UN
dat <- dat |> 
  mutate(UN = case_when(Q2_8 >= 4 ~ "Like",
                        Q2_8 == 3 ~ "Neither",
                        Q2_8 <= 2 ~ "Dislike"))
dat$UN <- as.factor(dat$UN)

# Shia party supporter
dat <- dat |> 
  mutate(Shia_party_supporter = (Q4_2 + Q4_3)/2) |>
  mutate(Shia_party_supporter = ifelse(Shia_party_supporter >= 7, "Shia-party supporter", "Others"))
dat$Shia_party_supporter <- as.factor(dat$Shia_party_supporter)

# Shia parties supporter Amal or Hezbollah: Alternative calculation way
#dat <- dat |> 
#  mutate(Shia_party_supporter = case_when(
#    Q4_2 >= 7 | Q4_3 >= 7 ~ "Shia Party Supporter",
    TRUE ~ "Others"
#  ))

#dat$Shia_party_supporter <- as.factor(dat$Shia_party_supporter)
#table(dat$Shia_party_supporter)

# Patronage
dat <- dat |> 
  mutate(Patronage = case_when(Q4 >= 3 ~ "High",
                               Q4 <= 2 ~ "Low"))
dat$Patronage <- as.factor(dat$Patronage)

# High patronage and Shia party supporter
dat <- dat |>
  mutate(Patron_Shia = ifelse(Shia_party_supporter == "Shia-party supporter" & Patronage == "High", "Patron_Shia", "Other"))
dat$Patron_Shia <- as.factor(dat$Patron_Shia)

# Rely job on state apparatuses
dat <- dat |> 
  mutate(Job_rely_state = ifelse(Q5_1 <= 5, "Rely job on state", "Others"))
dat$Job_rely_state <- as.factor(dat$Job_rely_state)

# Rely job on Zaim
dat <- dat |> 
  mutate(Job_rely_Zaim = ifelse(Q5_1 == 6, "Rely job on Zaim", "Others"))
dat$Job_rely_Zaim <- as.factor(dat$Job_rely_Zaim)

# Rely economic difficulty on state apparatuses
dat <- dat |> 
  mutate(economy_rely_state = ifelse(Q5_1 <= 5, "Rely economy on state", "Others"))
dat$economy_rely_state <- as.factor(dat$economy_rely_state)

# Rely security on state apparatuses
dat <- dat |> 
  mutate(Security_rely_state = ifelse(Q5_3 <= 5, "Rely security on state", "Others"))
dat$Security_rely_state <- as.factor(dat$Security_rely_state)

# Piety
dat <- dat |> 
  mutate(Piety = case_when(Q6 <= 4 ~ "Low",
                           Q6 >= 5 & Q6 <= 6 ~ "Middle",
                           Q6 >=7 ~ "High"))
dat$Piety <- as.factor(dat$Piety)

# Subjective conflict intensity
dat <- dat |> 
  mutate(Conflict_assessment = case_when(Q14 >= 4 ~ "Worse",
                                         Q14 == 3 ~ "Same",
                                         Q14 <= 2 ~ "Better"))
dat$Conflict_assessment <- as.factor(dat$Conflict_assessment)

# Subjective security perspective
dat$Q15 <- as.numeric(dat$Q15)
dat <- dat |> 
  mutate(Future_expectation = case_when(Q15 >= 4 ~ "Worse",
                                        Q15 == 3 ~ "Same",
                                        Q15 <= 2 ~ "Better"))
dat$Future_expectation <- as.factor(dat$Future_expectation)

```

# AMCE whole model
```{r, fig.height = 5, dpi = 300}
f1 <- selected ~ Activity + Armament + Patron + Service
str(dat)

# MM
whole_mm <- cj(dat, f1, 
               id = id,
               estimate = "mm", 
               h0 = 0.5,
               alpha = 0.1)
whole_mm
plot(whole_mm, vline = 0.5, 
     size = 2)+
  ggplot2::labs(title = "Whole result of MM")

whole_Lebanon <- plot(whole_mm, vline = 0, size = 2)+
  ggplot2::labs(title = "Whole result of AMCE")
ggsave(filename = "result/whole_MM.png",
       plot = whole_Lebanon,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

write_html(whole_mm, file = "result/whole_mm_Lebanon.html")

# amce
whole_amce <- cj(data = dat, f1,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 alpha = 0.1)
whole_amce
plot(whole_amce, vline = 0,
     size = 2) +
  ggplot2::labs(title = "Whole result of AMCE")

whole_Lebanon <- plot(whole_amce, vline = 0, size = 2)+
  ggplot2::labs(title = "Whole result of AMCE")
ggsave(filename = "result/whole_AMCE.png",
       plot = whole_Lebanon,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

write_html(whole_amce, file = "result/whole_amce_Lebanon.html")

# plot estimation difference
whole_mm$Estimate <- "MM"
whole_amce$Estimate <- "AMCE"
plot(rbind(whole_mm, whole_amce), feature_headers = FALSE, size = 2)+
  ggplot2::facet_wrap(~ Estimate, ncol = 2L)
```

# AMCE whole Model refined BW plot
```{r, fig.height = 5, dpi = 300}
shapes <- c("Activity" = 16, "Armament" = 17, "Patron" = 15, "Service" = 18)
whole_mm <- cj(dat, f1, 
               id = ~ respondentIndex, 
               estimate = "mm", 
               h0 = 0.5,
               alpha = 0.1)

plot_whole_mm <- plot(whole_mm, vline = 0.5, size = 2) +
  ggplot2::aes(shape = feature) +
  ggplot2::scale_shape_manual(values = shapes) +
  ggplot2::scale_color_grey(start = 0, end = 0) +
  ggplot2::theme_bw() +
  ggplot2::labs(title = "Whole result of MM", shape = "Attribute")

ggsave(filename = "result/whole_MM_bw.png",
       plot = plot_whole_mm, width = 7, height = 5, dpi = 600)

whole_amce <- cj(data = dat, f1,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 alpha = 0.1)

plot_whole_amce <- plot(whole_amce, vline = 0, size = 2) +
  ggplot2::aes(shape = feature) +
  ggplot2::scale_shape_manual(name = "Attribute", values = shapes) + 
  ggplot2::scale_color_grey(start = 0, end = 0) +
  ggplot2::theme_bw() +
  ggplot2::guides(
    color = "none", 
    shape = guide_legend(title = "Attribute")
  ) +
  ggplot2::labs(title = "Whole result of AMCE")

ggsave(filename = "result/whole_AMCE_bw.png",
       plot = plot_whole_amce, width = 7, height = 5, dpi = 600)

whole_mm$Estimate <- "MM"
whole_amce$Estimate <- "AMCE"

plot_comparison <- plot(rbind(whole_mm, whole_amce), feature_headers = FALSE, size = 2) +
  ggplot2::facet_wrap(~ Estimate, ncol = 2L) +
  ggplot2::aes(shape = feature) +
  ggplot2::scale_shape_manual(values = shapes) +
  ggplot2::scale_color_grey(start = 0, end = 0) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "bottom")

print(plot_comparison)
```

## H1a/b
# Shia party supporter
```{r, fig.height = 5, dpi = 300}
# 1. AMCE
Shia_party_supporter_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Shia_party_supporter,
                 alpha = 0.1)
Shia_party_supporter_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Shia_party_supporter,
                      alpha = 0.1)
plot(rbind(Shia_party_supporter_amce, Shia_party_supporter_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# 2. MM
Shia_party_supporter_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Shia_party_supporter,
                 alpha = 0.1)
Shia_party_supporter_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Shia_party_supporter,
                      alpha = 0.1)
plot(rbind(Shia_party_supporter_mm, Shia_party_supporter_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Shia_party_supporter_mm, group = "Shia_party_supporter", size = 2)


# MM plot
plot(Shia_party_supporter_mm, group = "Shia_party_supporter", size = 2)
Shia_party_supporter <- plot(Shia_party_supporter_mm, group = "Shia_party_supporter", size = 2) +
  ggplot2::labs(title = "Shia-party supporter")
ggsave(filename = "result/Shia_party_supporter_mm.png",
       plot = Shia_party_supporter,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

# MM difference plot
plot(Shia_party_supporter_mm, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

Shia_party_supporter_diff <- plot(Shia_party_supporter_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L) +
  ggplot2::labs(title = "Shia-party supporter")
ggsave(filename = "result/Shia_party_supporter_mm_diff.png",
       plot = Shia_party_supporter_diff,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

# anova
cj_anova(data = dat, 
         selected ~ Activity + Armament + Patron + Service, by = ~ Shia_party_supporter)

# 3. OLS result
write_html(Shia_party_supporter_mm, file = "result/Shia_party_supporter_mm.html")
write_html(Shia_party_supporter_mm_diff, file = "result/Shia_party_supporter_mm_diff.html")
```

# Shia party supporter refined BW plot
```{r, fig.height = 5, dpi = 300}
shapes <- c("Activity" = 16, "Armament" = 17, "Patron" = 15, "Service" = 18)

dat$Shia_party_supporter <- as.factor(dat$Shia_party_supporter)

Shia_party_supporter_amce <- cj(data = dat,
                                selected ~ Activity + Armament + Patron + Service,
                                id = ~ respondentIndex,
                                estimate = "amce",
                                by = ~ Shia_party_supporter,
                                alpha = 0.1)

Shia_party_supporter_amce_diff <- cj(data = dat,
                                     selected ~ Activity + Armament + Patron + Service,
                                     id = ~ respondentIndex,
                                     estimate = "amce_diff",
                                     by = ~ Shia_party_supporter,
                                     alpha = 0.1)

plot(rbind(Shia_party_supporter_amce, Shia_party_supporter_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L) +
  ggplot2::aes(shape = feature) +
  ggplot2::scale_shape_manual(name = "Attribute", values = shapes) +
  ggplot2::scale_color_manual(values = c("1" = "black", "0" = "gray60", 
                                         "TRUE" = "black", "FALSE" = "gray60")) +
  ggplot2::theme_bw() +
  ggplot2::guides(color = "none")

Shia_party_supporter <- plot(Shia_party_supporter_mm, group = "Shia_party_supporter", size = 2) +
  ggplot2::aes(shape = feature) +
  ggplot2::scale_shape_manual(values = shapes, name = "Attribute") +
  ggplot2::scale_color_manual(
    values = c("Shia-party supporter" = "black", "Shia Supporter" = "black", "1" = "black",
               "Others" = "gray60", "0" = "gray60"),
    na.translate = FALSE
  ) +
  
  ggplot2::theme_bw() +
  ggplot2::labs(title = "Shia-party supporter")

ggsave(filename = "result/Shia_party_supporter_mm.png",
       plot = Shia_party_supporter,
       width = 7, height = 4, units = "in", dpi = 600)

Shia_party_supporter_diff <- plot(Shia_party_supporter_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L) +
  ggplot2::aes(shape = feature) +
  ggplot2::scale_shape_manual(values = shapes, name = "Attribute") +
  ggplot2::scale_color_grey(start = 0, end = 0) + 
  ggplot2::theme_bw() +
  ggplot2::guides(color = "none") +
  ggplot2::labs(title = "Shia-party supporter (Difference)")

ggsave(filename = "result/Shia_party_supporter_mm_diff.png",
       plot = Shia_party_supporter_diff,
       width = 7, height = 4, units = "in", dpi = 600)
```

# High patronage and Shia_party_supporter
```{r, fig.height = 5, dpi = 300}
# amce
Patron_Shia_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Patron_Shia,
                 alpha = 0.1)
Patron_Shia_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Patron_Shia,
                      alpha = 0.1)
plot(rbind(Patron_Shia_amce, Patron_Shia_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Patron_Shia_mm <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Patron_Shia,
                 alpha = 0.1)
Patron_Shia_mm_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Patron_Shia,
                      alpha = 0.1)
plot(rbind(Patron_Shia_mm, Patron_Shia_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Patron_Shia_mm, group = "Patron_Shia", size = 2)

# MM plot
plot(Patron_Shia_mm, group = "Patron_Shia", size = 2)
Patron_Shia <- plot(Patron_Shia_mm, group = "Patron_Shia", size = 2) +
  ggplot2::labs(title = "Shia-party Patronage")
ggsave(filename = "result/Patron_Shia_mm.png",
       plot = Patron_Shia,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

# MM difference plot
plot(Patron_Shia_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

Patron_Shia_diff_LB <- plot(Patron_Shia_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L) +
  ggplot2::labs(title = "Shia-party Patronage")
ggsave(filename = "result/Patron_Shia_mm_diff.png",
       plot = Patron_Shia_diff_LB,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

# anova
cj_anova(data = dat, 
         selected ~ Activity + Armament + Patron + Service, by = ~ Patron_Shia)

# OLS result
write_html(Patron_Shia_mm, file = "result/Patron_Shia_mm.html")
write_html(Patron_Shia_mm_diff, file = "result/Patron_Shia_mm_diff.html")
```

# High patronage and Shia_party_supporter refined BW plot
```{r, fig.height = 5, dpi = 300}
shapes <- c("Activity" = 16, "Armament" = 17, "Patron" = 15, "Service" = 18)

dat$Patron_Shia <- factor(dat$Patron_Shia, 
                          levels = c("Patron_Shia", "Other"), 
                          labels = c("High", "Low"))

Patron_Shia_amce <- cj(data = dat,
                       selected ~ Activity + Armament + Patron + Service,
                       id = ~ respondentIndex,
                       estimate = "amce",
                       by = ~ Patron_Shia,
                       alpha = 0.1)

Patron_Shia_amce_diff <- cj(data = dat,
                            selected ~ Activity + Armament + Patron + Service,
                            id = ~ respondentIndex,
                            estimate = "amce_diff",
                            by = ~ Patron_Shia,
                            alpha = 0.1)

plot(rbind(Patron_Shia_amce, Patron_Shia_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L) +
  ggplot2::aes(shape = feature) +
  ggplot2::scale_shape_manual(name = "Attribute", values = shapes) +
  ggplot2::scale_color_manual(values = c("High" = "black", "Low" = "gray60")) +
  ggplot2::theme_bw() +
  ggplot2::guides(color = "none")

Patron_Shia_mm <- cj(data = dat,
                     selected ~ Activity + Armament + Patron + Service,
                     id = ~ respondentIndex,
                     estimate = "mm",
                     by = ~ Patron_Shia,
                     alpha = 0.1)

Patron_Shia_mm_diff <- cj(data = dat,
                          selected ~ Activity + Armament + Patron + Service,
                          id = ~ respondentIndex,
                          estimate = "mm_diff",
                          by = ~ Patron_Shia,
                          alpha = 0.1)

Patron_Shia_plot <- plot(Patron_Shia_mm, group = "BY", size = 2) +
  ggplot2::aes(shape = feature) +
  ggplot2::scale_shape_manual(name = "Attribute", values = shapes) +
  ggplot2::scale_color_manual(
    name = "Group",
    values = c("High" = "black", "Low" = "gray60"),
    na.translate = FALSE
  ) + 
  
  ggplot2::theme_bw() +
  ggplot2::labs(title = "Shia-party Patronage") +
  
  ggplot2::guides(
    fill = "none",
    color = guide_legend(title = "Group"),
    shape = guide_legend(title = "Attribute")
  ) +
  ggplot2::theme(
    legend.position = "right",
    legend.box.margin = margin(l = 5)
  )

ggsave(filename = "result/Patron_Shia_mm.png",
       plot = Patron_Shia_plot,
       width = 8, height = 5, units = "in", dpi = 600)

Patron_Shia_diff_LB <- plot(Patron_Shia_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L) +
  ggplot2::aes(shape = feature) +
  ggplot2::scale_shape_manual(name = "Attribute", values = shapes) +
  ggplot2::scale_color_grey(start = 0, end = 0) + 
  ggplot2::theme_bw() +
  ggplot2::guides(color = "none") +
  ggplot2::labs(title = "Shia-party Patronage (Difference)") +
  
  ggplot2::theme(
    legend.position = "right",
    legend.box.margin = margin(l = 5)
  )

ggsave(filename = "result/Patron_Shia_mm_diff.png",
       plot = Patron_Shia_diff_LB,
       width = 8, height = 5, units = "in", dpi = 600)

cj_anova(data = dat, 
         selected ~ Activity + Armament + Patron + Service, by = ~ Patron_Shia)

write_html(Patron_Shia_mm, file = "result/Patron_Shia_mm.html")
write_html(Patron_Shia_mm_diff, file = "result/Patron_Shia_mm_diff.html")
```

## H2
# Ceasefire
```{r, fig.height = 5, dpi = 300}
# 1. AMCE
ceasefire_amce <- cj(data = dat,
                 selected ~ Activity + Armament + Patron + Service,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ ceasefire,
                 alpha = 0.1)
ceasefire_amce_diff <- cj(data = dat,
                      selected ~ Activity + Armament + Patron + Service,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ ceasefire,
                      alpha = 0.1)
plot(rbind(ceasefire_amce, ceasefire_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# 2. MM
dat$ceasefire <- relevel(dat$ceasefire, "Pre")
ceasefire_mm <- cj(dat,
               selected ~ Activity + Armament + Patron + Service,
               id = id,
               estimate = "mm",
               by = ~ ceasefire,
               alpha = 0.1)
ceasefire_mm_diff <- cj(data = dat,
                    selected ~ Activity + Armament + Patron + Service,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ ceasefire,
                    alpha = 0.1)
plot(rbind(ceasefire_mm, ceasefire_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)


# MM plot
plot(ceasefire_mm, group = "ceasefire", size = 2)
ceasefire <- plot(ceasefire_mm, group = "ceasefire", size = 2) +
  ggplot2::labs(title = "Ceasefire")
ggsave(filename = "result/ceasefire.mm.png",
       plot = ceasefire,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

# MM difference plot
plot(ceasefire_mm, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

ceasefire_diff <- plot(ceasefire_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L) +
  ggplot2::labs(title = "Ceasefire")
ggsave(filename = "result/ceasefire_mm_diff.png",
       plot = ceasefire_diff,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

# 3. OLS
write_html(ceasefire_mm, file = "result/ceasefire_mm.html")
write_html(ceasefire_mm_diff, file = "result/ceasefire_mm_diff.html")

# anova_nested model
cj_anova(data = dat, 
         selected ~ Activity + Armament + Patron + Service, by = ~ ceasefire)
```

# Ceasefire refined BW plot
```{r, fig.height = 5, dpi = 300}
shapes <- c("Activity" = 16, "Armament" = 17, "Patron" = 15, "Service" = 18)

dat$ceasefire <- as.factor(dat$ceasefire)
dat$ceasefire <- relevel(dat$ceasefire, "Pre")

ceasefire_mm <- cj(data = dat,
                   selected ~ Activity + Armament + Patron + Service,
                   id = ~ respondentIndex,
                   estimate = "mm",
                   by = ~ ceasefire,
                   alpha = 0.1)

ceasefire_mm_diff <- cj(data = dat,
                        selected ~ Activity + Armament + Patron + Service,
                        id = ~ respondentIndex,
                        estimate = "mm_diff",
                        by = ~ ceasefire,
                        alpha = 0.1)

ceasefire_plot <- plot(ceasefire_mm, group = "BY", size = 2) +
  ggplot2::aes(shape = feature) +
  ggplot2::scale_shape_manual(name = "Attribute", values = shapes) +
  ggplot2::scale_color_grey(
    start = 0, end = 0.6, 
    na.translate = FALSE 
  ) + 
  
  ggplot2::theme_bw() +
  ggplot2::labs(title = "Ceasefire (Marginal Means)") +
  ggplot2::guides(
    fill = "none",
    color = guide_legend(title = "Ceasefire Status"),
    shape = guide_legend(title = "Attribute")
  ) +
  ggplot2::theme(
    legend.position = "right",
    legend.box.margin = margin(l = 5) 
  )

ggsave(filename = "result/ceasefire_mm.png",
       plot = ceasefire_plot,
       width = 8, height = 5, units = "in", dpi = 600)

ceasefire_diff_plot <- plot(ceasefire_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L) +
  ggplot2::aes(shape = feature) +
  ggplot2::scale_shape_manual(name = "Attribute", values = shapes) +
  ggplot2::scale_color_grey(start = 0, end = 0) + 
  ggplot2::theme_bw() +
  ggplot2::guides(color = "none") +
  ggplot2::labs(title = "Ceasefire (Difference)") +
  ggplot2::theme(
    legend.position = "right",
    legend.box.margin = margin(l = 5)
  )

ggsave(filename = "result/ceasefire_mm_diff.png",
       plot = ceasefire_diff_plot,
       width = 8, height = 5, units = "in", dpi = 600)
```

## Appendix additional analyses
# Interaction: ceasefire*Shia_party_supporter_DD
```{r, fig.height = 5, dpi = 300}
dat$Shia_party_ceasefire <- interaction(dat$Shia_party_supporter, dat$ceasefire, sep = "_")
Shia_party_ceasefire_mm <- cj(dat,
               selected ~ Shia_party_ceasefire,
               id = id,
               estimate = "mm",
               alpha = 0.1)
plot(Shia_party_ceasefire_mm, vline = 0.5, size = 2)

# plot
Shia_party_ceasefire <- plot(Shia_party_ceasefire_mm, vline = 0, size = 2)+
  labs(title = "Shia-party supporter x ceasefire")
ggsave(filename = "result/Shia_party_ceasefire_mm.png",
       plot = Shia_party_ceasefire,
       width = 5,
       height = 3,
       units = "in", 
       dpi = 600)

# OLS result
write_html(Shia_party_ceasefire_mm, file = "result/Shia_party_ceasefire_mm.html")
```

# Interaction: ceasefire*Shia_party_supporter_DD refined
```{r, fig.height = 5, dpi = 300}
dat$ceasefire <- factor(dat$ceasefire, levels = c("Pre", "Post"))
dat$Shia_party_ceasefire <- interaction(dat$Shia_party_supporter, dat$ceasefire, sep = "_")

Shia_party_ceasefire_mm <- cj(dat,
                              selected ~ Shia_party_ceasefire,
                              id = ~ respondentIndex,
                              estimate = "mm",
                              alpha = 0.1)

plot_data <- Shia_party_ceasefire_mm %>%
  mutate(
    Group = case_when(
      grepl("Shia-party supporter", level) ~ "Shia-party Supporter",
      grepl("Others", level) ~ "Others"
    ),
    Period = case_when(
      grepl("Pre", level) ~ "Pre-Ceasefire",
      grepl("Post", level) ~ "Post-Ceasefire"
    )
  ) %>%
  mutate(level = factor(level, levels = c(
    "Others_Pre", "Others_Post", 
    "Shia-party supporter_Pre", "Shia-party supporter_Post"
  )))

p_interaction <- ggplot(plot_data, aes(x = estimate, y = level, color = Group)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2, size = 0.8) +
  geom_point(size = 4) +
  scale_y_discrete(labels = c(
    "Others_Pre" = "Others (Pre)",
    "Others_Post" = "Others (Post)",
    "Shia-party supporter_Pre" = "Shia-party Supporter (Pre)",
    "Shia-party supporter_Post" = "Shia-party Supporter (Post)"
  )) +
  scale_color_manual(values = c("Others" = "gray60", "Shia-party Supporter" = "black")) +
  labs(title = "Support Level by Group & Period",
       subtitle = "Interaction: Supporter Status × Ceasefire",
       x = "Marginal Mean (Support)",
       y = "") +
  theme_bw(base_family = "sans", base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    axis.text.y = element_text(color = "black", size = 12),
    axis.text.x = element_text(color = "black"),
    legend.position = "none", 
    panel.grid.major.y = element_blank(), 
    panel.grid.major.x = element_line(color = "gray90", linetype = "dotted")
  )

print(p_interaction)

ggsave(filename = "result/Shia_party_ceasefire_mm.png",
       plot = p_interaction,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

write_html(Shia_party_ceasefire_mm, file = "result/Shia_party_ceasefire_mm.html")
```

# Interaction: Shia party supporter x ceasefire
```{r}
table(dat$ceasefire, dat$Shia_party_supporter)
cor(as.numeric(dat$ceasefire), as.numeric(dat$Shia_party_supporter))

# subset ceasefire and shia party supporter
pre_supporter <- dat |>
  filter(ceasefire == "Pre")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ Shia_party_supporter,
     alpha = 0.1)
plot(pre_supporter, group = "Shia_party_supporter", 
     size = 2,
     vline = 0.5)

post_supporter <- dat |>
  filter(ceasefire == "Post")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ Shia_party_supporter,
     alpha = 0.1)
plot(post_supporter, group = "Shia_party_supporter", 
     size = 2,
     vline = 0.5)

plot(rbind(pre_supporter, post_supporter), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

pre_supporter_diff <- dat |>
  filter(ceasefire == "Pre")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ Shia_party_supporter,
     alpha = 0.1)
plot(pre_supporter_diff, group = "Shia_party_supporter", 
     size = 2,
     vline = 0.5)

post_supporter_diff <- dat |>
  filter(ceasefire == "Post")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ Shia_party_supporter,
     alpha = 0.1)
plot(post_supporter_diff, group = "Shia_party_supporter", 
     size = 2,
     vline = 0.5)

plot(rbind(pre_supporter, pre_supporter_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
plot(rbind(post_supporter, post_supporter_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)


plot(rbind(pre_supporter_diff, post_supporter_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)


###
dat$ceasefire <- relevel(dat$ceasefire, "Pre")
supporter <- dat |>
  filter(Shia_party_supporter == "Shia-party supporter")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(supporter, group = "ceasefire", 
     size = 2,
     vline = 0.5)

non_supporter <- dat |>
  filter(Shia_party_supporter == "Others")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(non_supporter, group = "ceasefire", 
     size = 2,
     vline = 0.5)
plot(rbind(supporter, non_supporter), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

supporter_diff <- dat |>
  filter(Shia_party_supporter == "Shia-party supporter")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(supporter_diff, group = "ceasefire", 
     size = 2)

non_supporter_diff <- dat |>
  filter(Shia_party_supporter == "Others")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(non_supporter_diff, group = "ceasefire", 
     size = 2)

# interaction
supporter$Supporter <- "Supporter"
non_supporter$Supporter <- "Non_supporter"
interaction <- rbind(supporter, non_supporter)

H1_app1 <- plot(interaction, group = "Supporter", size = 2) +
  ggplot2::facet_wrap(~ BY, ncol = 2) +
  labs(title = "Shia-party supporter",
       shape = "Hezbollah Supporter",
       labels = c("Non-suporter", "Supporter"))+
  scale_color_hue(labels = c("Non-suporter", "Supporter"))
H1_app1

ggsave(filename = "result/H1_app1.png",
       plot = H1_app1,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)  

# OLS result
write_html(interaction, file = "result/H1_app1.html")

# diff
supporter_diff$Supporter <- "Supporter"
non_supporter_diff$Supporter <- "Non_supporter"
interaction_diff <- rbind(supporter_diff, non_supporter_diff)

H1_app2 <- plot(interaction_diff, group = "Supporter", size = 2) +
  ggplot2::facet_wrap(~ BY, ncol = 2) +
  labs(title = "Shia-party supporter",
       labels = c("Non-suporter", "Supporter"))+
  scale_color_hue(labels = c("Non-suporter", "Supporter"))
H1_app2

ggsave(filename = "result/H1_app2.png",
       plot = H1_app2,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)  
```

# additional
```{r}
table(dat$ceasefire, dat$Shia_party_supporter)
cor(as.numeric(dat$ceasefire), as.numeric(dat$Shia_party_supporter))

# subset ceasefire and shia party supporter
# mm
pre_supporter <- dat |>
  filter(ceasefire == "Pre")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ Shia_party_supporter,
     alpha = 0.1)
plot(pre_supporter, group = "Shia_party_supporter", 
     size = 2,
     vline = 0.5)

post_supporter <- dat |>
  filter(ceasefire == "Post")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ Shia_party_supporter,
     alpha = 0.1)
plot(post_supporter, group = "Shia_party_supporter", 
     size = 2,
     vline = 0.5)

plot(rbind(pre_supporter, post_supporter), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)


# mm diff
pre_supporter_diff <- dat |>
  filter(ceasefire == "Pre")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ Shia_party_supporter,
     alpha = 0.1)
plot(pre_supporter_diff, group = "Shia_party_supporter", 
     size = 2,
     vline = 0.5)

post_supporter_diff <- dat |>
  filter(ceasefire == "Post")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ Shia_party_supporter,
     alpha = 0.1)
plot(post_supporter_diff, group = "Shia_party_supporter", 
     size = 2,
     vline = 0.5)

plot(rbind(pre_supporter, pre_supporter_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
plot(rbind(post_supporter, post_supporter_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)


plot(rbind(pre_supporter_diff, post_supporter_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)


# mm+mm diff
dat$ceasefire <- relevel(dat$ceasefire, "Pre")
supporter <- dat |>
  filter(Shia_party_supporter == "Shia-party supporter")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(supporter, group = "ceasefire", 
     size = 2,
     vline = 0.5)

non_supporter <- dat |>
  filter(Shia_party_supporter == "Others")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(non_supporter, group = "ceasefire", 
     size = 2,
     vline = 0.5)
plot(rbind(supporter, non_supporter), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

supporter_diff <- dat |>
  filter(Shia_party_supporter == "Shia-party supporter")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(supporter_diff, group = "ceasefire", 
     size = 2)

non_supporter_diff <- dat |>
  filter(Shia_party_supporter == "Others")  |>
  cj(f1, 
     id = ~ respondentIndex, 
     estimate = "mm_diff", 
     by = ~ ceasefire,
     alpha = 0.1)
plot(non_supporter_diff, group = "ceasefire", 
     size = 2)

# OLS result
write_html(pre_supporter, post_supporter, file = "result/pre_post_supporter_mm.html")
write_html(pre_supporter_diff, post_supporter_diff, file = "result/pre_post_supporter_mm_diff.html")
```
